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Abstract 


The mechanical behavior of materials is usually simulated by the continuous mechanics approach. However, simula- 
tion of non-continuous phenomena like multi-fracturing is not well adapted to a continuous description. In this case, 
the discrete element method (DEM) is a good alternative because it naturally takes into account discontinuities. 
Many researchers have shown interest in this approach for wear and fracture simulation. The problem is that, while 
DEM is well adapted to simulate discontinuities, it is not suitable to simulate continuous behavior. In problems of 
wear or fracture, material is composed of continuous parts and discontinuous interfaces. The aim of the present work 
is to improve the ability of DEM to simulate the continuous part of the material using cohesive bond model. 
Continuous mechanics laws cannot be used directly within a DEM formulation. A second difficulty is that the 
volume between the discrete elements creates an artificial void inside the material. This paper proposes a methodology 
that tackles these theoretical difficulties and simulates, using a discrete element model, any material defined by a 
Young’s modulus, Poisson’s ratio and density, to fit the static and dynamic mechanical behavior of the material. The 
chosen cohesive beam model is shown to be robust concerning the influence of the discrete element sizes. This method 


is applied to a material which can be considered as perfectly elastic: fused silica. 


Keywords: Discrete element method, DEM, Calibration, Elastic, Dynamic, Fused silica 


1. Introduction 


The discrete element method (DEM) can describe 
quite naturally a granular medium. However, the num- 
ber of discrete elements to manage is high and it re- 
quires computational resources. The development of 
this method began in the early 1980s (1). More recently, 
researchers have used this method to study the damage 
of heterogeneous solids such as concrete (2) or rock (3), 
and homogeneous materials such as ceramics (4). 

Discrete element model is well adapted to simulate 
a media that has a great number of interfaces. It has 
been widely used to study tribological problems like 
wear phenomena. In this kind of problem, the mate- 
rial has a continuous part (the volume above the sur- 
face that is not yet affected by the wear), a continuous 
part with cracks (called subsurface damage in abrasion 
process terminology) and a discontinuous part (the in- 
terfacial media, called third body, that is a mixture of 
abrasive particles and wear particles). Discrete element 
model must be able to simulate with accuracy all these 
parts of the material. Unlike continuous approaches, the 
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main difficulty for DEM is to simulate properly the con- 
tinuous material. 

This paper focuses on a material which can be con- 
sideblue as homogeneous, isotropic and perfectly elas- 
tic: fused silica. This work is a continuation of a previ- 
ous study that investigated subsurface damage in silica 
glass due to surface polishing (5). In this previous study, 
discrete element models have shown qualitatively good 
agreement with experiments. The challenge, now, is to 
propose a 3D DEM spherical model allowing quanti- 
tative results for the silica considered as a continuous 
media. 

A preliminary task is the cohesive bond choice to 
model correctly the subsurface damage layer during the 
silica abrasive process. There are two main approaches 
to DEM : the dual spring model (a pair of normal and 
tangential springs) (6, 7, 8, 9) and the cohesive beam 
model (10, 11, 12). The beam cohesive model is not 
as well established in the literature as the classical dual 
spring model. Therefore, Schlangen and Garboczi in 
(10, §3) have shown that the beam cohesive bond pro- 
duces more realistic crack pattern than the simple spring 
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model or the dual spring model . A further work of these 
authors (11, $3) have studied the influence of the ini- 
tial geometric arrangement. They conclude that a disor- 
dered configuration gives more realistic crack geometry 
than an ordered geometrical configuration. 

According to these results, the discrete element ap- 
proach used in this paper is quite similar to the model 
described by Carmona in (12). An initial spherical do- 
main is generated by a numerical compaction method. 
A dispersion applied to the discrete element radii allows 
random geometrical arrangement. Then, the discrete el- 
ements are connected by the "beam cohesive bonds". 

The proposed model must simulate static and dy- 
namic behavior characterized by Young’s modulus, 
Poisson’s ratio, mechanical wave celerity and natural 
frequencies. 

The difficulty is that the mechanical behavior of a 
structure composed of a large number of discrete el- 
ements cannot be analytically predicted. Global be- 
haviors are the result of a large number of interactions 
between discrete elements and can be considered as 
an emergent physical property (13, preface). Implic- 
itly, two scales are considered in a discrete element ap- 
proach: 


e the structure scale, represented by a set of discrete 
elements. This scale will be called the "macro- 
scopic scale". 


e the discrete element scale and its elementary inter- 
action with its neighbors. This scale will be called 
the "microscopic scale". 


Note that the interesting properties, e.g., the Young’s 
modulus, the Poisson’s ratio, mechanical wave celer- 
ity and natural frequencies should be considered, in the 
DEM model, as emergent properties at the macroscopic 
scale. Furthermore, unlike the finite element method 
(14), continuous mechanical behavior laws cannot be 
directly introduced into the DEM formulation. As a re- 
sult, the difficulty is to quantify DEM microscopic in- 
teraction laws according to continuous mechanical be- 
havior. 

This problema has been discussed in detail by Ostoja- 
Starzewski (15). The author proposes micro-macro laws 
for some typical ordered lattice configurations. In the 
last paragraph (§ 6.3) dedicated to the periodic random 
lattice network, Ostoja-Starzewski proposes numerical 
tests to calibrate the model. The analytical approach is 
limited to an ordered and homogeneous configuration. 
This idea is well synthesized by Potyondy and Cundall 
(6, §3.1) who write: 

"For continuum models, the input properties (such as 


modulus and strength) can be derived directly from 
measurements performed on laboratory specimens. For 
the BPM! |...) the input properties of the components 
usually are not known. (...) For the general case of 
arbitrary packing of arbitrarily sized particles, the re- 
lation is found by means of a calibration process (...)". 
To summarize, for the random domains, the quantifica- 
tion of the microscopic parameters requires some nu- 
merical tests called a calibration procedure. This prob- 
lem has been intensively studied for the cohesive dual 
spring model. Hentz et al. (7) have proposed numer- 
ical quasi-static uni-axial traction tests to calibrate the 
bond parameters in regard to the macroscopic Young’s 
modulus and Poisson’s ratio. Note that these authors 
have introduced an energy criterion to reduce dispersion 
of the macroscopic emergent properties. Then, numer- 
ical dynamic tests are used to check the dynamic prop- 
erties. Fakhimi et al. (8) have used calibration curves 
and dimensional analysis instead of the trial and error 
approach. Tamarez and Plesha (9) have used analytical 
formulations of a 2D elementary volume (a unit cell). 

The cohesive beam model was first introduced by 
Herrmann in 1988 (16). This model was used in a 2D 
ordered lattice network (17, 18), then in a disordered 2D 
lattice network (11, 19, 20, 21, 22). In reference (11), 
the authors have considered that the microscopic and 
macroscopic Young’s moduli and Poisson’s ratios must 
be similar. The Beam dimensions (cross section and in- 
ertia momentum) are chosen thanks to a numerical re- 
cursive algorithm to satisfy an uniform elastic contin- 
uum condition. The main subject of other papers is the 
development of fracture models. The calibration meth- 
ods are not described in depth. Researchers simply rec- 
ommend using experimental and numerical tests. 

In this study, unlike that of Schlangen(11), the me- 
chanical properties of the cohesive beams will not be 
considered as similar to the reference material. The mi- 
croscopic local properties could be driven to produce 
discrete matter internal reorganization under loading 
that induces the desired behaviors at the macroscopic 
scale. 

In this paper the 3D cohesive beam model is briefly 
described in the first section. Next, a new calibration 
method adapted to the simulation of a perfectly homo- 
geneous, isotropic, elastic material with 3D discrete el- 
ements bonded by beams is introduced in three steps. 


1. The first step is the geometrical analysis of the ini- 
tial domain configuration, presented in three parts. 


!"Bonded Particle Model", discrete element model used by Po- 
tyondy and Cundall 


homogeneity: Voids inside the discrete material 
are minimized thanks to a compaction algo- 
rithm. A method is developed to ensure the 
validity of the sphere packing with appropri- 
ate control criteria. 

isotropy: An original and simple method to define 
and quantify the geometrical isotropy is pre- 
sented. 


fineness: A solution to define the micro-scale and 
macro-scale ratio (the critical number of dis- 
crete elements number that allows the emer- 
gence of stable geometrical properties) is dis- 
cussed. 


2. Then the influence of the beam parameters on the 
macroscopic behaviors is studied by quasi-static 
simulations. A micro-macro tendency is used to 
develop an original calibration method for the elas- 
tic parameters (macroscopic Young’s modulus and 
Poisson’s ratio). 

3. The same methodology is used for the calibration 
of the mass properties, by dynamic simulations. 

4. The influence of the discrete element size on the 
calibration results is studied. The independence 
of the calibrated parameters with regard to the dis- 
crete element size is a property that appears to be 
essential when quantitative results are wanted. 


The overall method is applied to the studied mate- 
rial, the silica glass. The calibration method is validated 
through a large number of static and dynamic tests. 


2. Cohesive beam bond model 


Figure | shows two discrete elements bonded by a 
cohesive beam. The cylindrical geometry is chosen be- 
cause it’s dimensional description requires only two in- 
dependent parameters: a length L, and a radius TE The 
mechanical properties are also linked to the cohesive 
beams: a Young’s modulus E, and a Poisson’s ratio vy. 
These four geometric and mechanical parameters suf- 
fice to describe the cohesive beam. Note that the cohe- 
sive beams are mass-less; mass properties are assigned 
only to the discrete elements. 

For the sake of clarity figure 2 shows a configuration 
in which the discrete elements have been moved away. 
The cohesive beam is symbolized by its median line. 
Both cohesive bond ends are fixed to the discrete ele- 
ment centers O; and O2. The discrete element frames 


?To distinguish micro from macro properties, micro parameters are 
denoted by °x’ and macro parameters by ° M’. 


Figure 1: The cohesive beam bond 
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Fi (01.%1.71.Z1) and F> (02,2, ¥2,2) are oriented 


such that x and x are normal to the beam cross section 
ends. At the initial time, the beams are relaxed (figure 
2a). Figure 2b shows the cohesive beam in a loading 
state induced by the discrete element movement rela- 
tively to the initial configuration. 


The analytical model of Euler-Bernoulli beam is well 
known (23). In reference (24, $6.2), the author describe 
a stiffness matrix expressed in the beam local frame for 
a finite element application. Figure 2b illustrates the 
beam local frame positioning. The center of discrete el- 
ement 1 (Q;) is considered as the origin. The "aligned" 

——> => > 
configuration, in which 0102 = kX; = —kX2, is con- 
sidered as the non-bending state and is taken as ref- 
erence. Consequently, the cohesive beam local frame 


>> CO). . 
F Q X,Y, 2) is oriented such that (see figure 2b): 


In the local frame F, the deflections at O; and O2 are 


null. Cross section bending rotations at O; and O2 
i> —> 


are defined, respectively, by 0; = (x Xi) and b = 


— > | 
(=x i X). Consequently, the force and torque reactions 


acting on discrete elements 1 and 2 are: 


(a) Relaxing state 


Discrete Element 1 


Discrete Element 2 


(b) Loading state 


Figure 2: Cohesive beam bond configuration 
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—> D ip N 
e Fpe1 and Fpp2 are the beam force reactions acting 
on discrete elements 1 and 2. 


——> — s 5 
e Tpg and T hp) are the beam torque reactions acting 
on discrete elements 1 and 2. 


e l, and A, are the initial beam length and the lon- 
gitudinal extension. 


= > : 
© A(A1x, Oiy, 12) and 42(62,, 2), 6.) are the rotations 


of beam cross section at the points O; and O; ex- 
pressed in the beam local frame. 


e S, 10, and 1, are the beam cross section area, po- 
lar moment of inertia and moment of inertia along 


=> => 
Y and Z. 
e E, and G, are the Young and shear modulus. 


Note that reaction force and torque are expressed in the 
= = = 
beam local frame F |O, X, Y, Z |. The four parameters 


that define the micro beam, L,, fu, Eu and v, have only 
an influence on the elastic behavior of the assembly. 


3. Explicit dynamic algorithm 


The numerical resolution is based on an explicit in- 
tegration scheme well adapted to massive DEM simu- 
lation (25) and high velocity phenomena such as frac- 
turing or impact simulation. Many explicit schemes can 
be used : the Verlet velocity, Runge-Kutta, leapfrog or 
gear’s method. .. (26, $13). In reference (25), the au- 
thors have compared these algorithms "in terms of ac- 
curacy, Stability and CPU efficiency". It appears that all 
of them give approximately the same efficiency. 

Velocity Verlet scheme is chosen for its simplicity . 
Discrete element position and velocity are estimated by: 


>, >, > dt > 

P(t + dt) = p(t) + dt pO + 7 PO (5) 
> > dt / > 2 

Be + an = PO + (PO + BE+dD) © 


where : 


e tis the current time and dt is the integration time 
step. 


e p(t), p(t) and p(t) is the discrete element linear po- 
sition, velocity and acceleration. 


The discrete element orientations are described by 
quaternions, noted q(t), that allow an efficiency way to 
compute the rotation of the local frames associated with 
discrete elements (27, §2.5). Quaternion is linked to the 
angular velocity with the following equality (26, $10.5): 


1 
q(t) = 700) q(t) (7) 


Where &(?) is the angular velocity of a discrete element. 
The velocity Verlet scheme is also applied to quaternion 
q(t), with: 


d? 
q(t + dt) = q(t) + dt q(t) + iO (8) 


gasd =d+ Oria O 


To prevent q(t) numerical drift, the quaternion must be 
normalized at each time step. Algorithm 1 details the 
encapsulation of Verlet velocity in an explicit dynamic 
resolution. Note that this numerical scheme is not well- 
adapted to quasi-static simulation. Special care, de- 
scribed later in this paper, will be taken with this kind 
of test. 


Algorithm 1 Explicit dynamic resolution 


Require: (0) PO) ËO) q(0) 4(0) 4(0) 
t— 0 
for all iteration n do 
for all discrete element i do 

Pi(t + dt) + velocity Verlet scheme (eq. 5) 
fit + dt) — Sum of force acting on i 
p(t + dt) — Newton second law 
p(t + dt) + velocity Verlet scheme (eq. 6) 


qi(t + dt) — Velocity Verlet scheme (eq. 8) 
qi(t + dt) — Normalization 
7(t + dt) — Sum of torque acting on i 
git + dt) — Angular momentum law 
gi(t + dt) — Velocity Verlet scheme (eq. 9) 
end for 
t<t+dt 
end for 


4. Geometrical properties of the initial discrete do- 
main 


The initial discrete domain geometry must be in ac- 
cordance with the structural properties of the simulated 


material. In the case of the perfectly elastic solid, the 
main properties are homogeneity and isotropy. The 
geometrical discrete domain configuration impacts the 
mechanical behavior at the macroscopic scale. In ref. 
(11, 12) the authors have observed that ordered arrange- 
ments promote anisotropic phenoma such as preferred 
crack paths or the non-uniform propagation of elastic 
waves. 

In the case of the spherical discrete elements, artifi- 
cial voids are generated between discrete elements in- 
side the material. At this scale, the discrete material is 
not homogeneous. To reduce undesirable voids, the dis- 
crete domain must be initially compacted. 

This subsection first discusses a method that ensures 
and verifies the geometrical homogeneity and isotropy. 
Then, the influence of the discrete domain refining is 
studied. 


4.1. Homogeneity 


To decrease the artificial voids between the discrete 
elements, various packing processes can be used, for 
example, dynamic methods such as iterative growth al- 
gorithm (28) or isotropic compression (29) or geometric 
algorithm (30). In this study, the discrete domain is gen- 
erated following a dynamic custom recipe. The figure 3 
shows the initial configuration. The discrete elements 
are placed along a plane-parallel grid. an uniform dis- 
persion law, detailed later in this section, is applied to 
discrete element radii. The compaction simulation con- 
sist of : 


e applying an horizontal sinusoidal movement to the 
"Shear Wall" set. 


e applying a vertical pressure to "Pressure Wall" set. 


e confining the discrete domain in a given volume 
thanks to the "Repulsive Walls" which apply an 
elastic repulsion law. 


For an identical or small dispersion radius, the sphere 
packing should be similar to a "Random Close Packing" 
(RCP) as described in (31). A first step to validate the 
initial discrete domain is to check its conformity with 
RCP definition. The domain is considered as correctly 
compacted, if the sphere packing gives a volume frac- 
tion value around 0.636 (32) and a cardinal number? 
value around 6 (33). However, these two values do not 
ensure that the packing is geometrically isotropic. 


3 Average number of contact per discrete element. 


Pressure 


Pressure Wall 


So 


Shear Wall displacement 


Figure 3: Shear compaction method 


4.2. Isotropy 


The definition of geometrical isotropy must be clar- 
ified before proposing a criterion. Cambou (34, Intro- 
duction, $3.6) defines the geometric anisotropy as the 
distribution of contact directions. If this distribution is 
perfectly homogeneous, the domain is considered as ge- 
ometrically isotropic. To summarize, an assertion could 
be formulated as: the geometrical isotropy is a neces- 
sary condition to ensure the mechanical isotropy of the 
simulated material. To "measure" the geometrical orga- 
nization of granular material, the authors in (35, $1.2.2) 
have exploited a mathematical tool called the "fabric 
tensor". However, this tool cannot be used to deter- 
mine the isotropy in a simple way (36). A more intu- 
itive method based on a simple geometric computation 
and statistical analysis is proposed. It is an extension 
to 3D space of classical 2D graphs that classify contact 
into a direction subset (for example references (36) and 
(20)). 

Contacts are grouped into a subset depending on their 
orientation in 3D space. All the members in a set have 
a quite similar spatial orientation. To group contacts 
in orientation sets, a platonic solid (a "geode") of 320 
equal faces is built. Each contact is placed at the geode 
center, and added to the corresponding intersected face 
group (see figure 4). The result is a kind of 3D his- 
togram in which each class maps an orientation subset 
(a sort of discretized solid angle); any class weight gives 
the density of contact orientations that matches the solid 
angle. For a low dispersion of class weights, the dis- 
crete domain is considered as geometrically isotropic 
(an homogeneous distribution of contact orientation in 
3D space). 


Figure 4: Platonic solid (Geode) used as a reference ge- 
ometry to classify the contact orientation 


To illustrate this method, the influence of the parti- 
cle size distribution on isotropy is studied. It is known 
from literature (31) that to prevent an ordered packing 
configuration, a dispersion (noted x) must be applied on 
discrete elements radius. In this study, the dispersion x 
is uniform and is defined as the dispersion range : 

_ Rinax = Rmin 


A (10) 


where Rmax, Rmin and R are the maximum, minimum and 
average discrete element radius values. 

Figures 5 shows the geometrical arrangement for two 
values of the dispersion parameter x. The figure illus- 
trates the influence of radius dispersion on the geomet- 
rical arrangement (figures 5a and 5b) and on the contact 
orientation (figures 5c and 5d). For a distribution range 
k = 0% the packing seems to be perfectly ordered. The 
perfect arrangement is strongly anisotropic. In contrast, 
a higher dispersion value (k = 25%) seems to promote 
the isotropy. 

Figures 6a and 6b show the 3D histograms used to 
qualify the observed level of isotropy. From these fig- 
ures, it is clear that the radius dispersion value x highly 
influences the isotropy level. To quantify this level, it 
is proposed to compute the mean square difference be- 
tween observed frequencies in geode cells (f;) and the 
uniform frequency (1/N) : 


N (pel à 
e= SET) x) (11) 


(b) Domain for xk = 25% 
(discrete element view) 


(a) Domain for x = 0% 
(discrete element view) 


(d) Domain for « = 25% 
(contact view) 


(c) Domain for x = 0% 
(contact view) 


Figure 5: Geometrical arrangement of a 3D sphere 
packing with different values of x 


Where : 


e N is the total number of cells, i.e, the 320 "geode" 
faces. 


e f; is the observed frequency of the cell i, i.e, the 
ratio between the number of contact that matches 
the solid angle i and the total contact number. 


The important aspect of this criterion is the asymptotic 
behavior (see figure 7). Increasing radius dispersion 
value x gives an asymptotic constant limit, beginning 
from a x value of 15 %. This result is in accordance 
with the observation of Luding (37, chapter 5): "crys- 
tallization (... ) does not occur for polydisperse packing 
with wo Z 0.15”. In other words, for a dispersion value 
higher than 15%, an ordered geometrical arrangement 
does not occur within the sphere packing. Thus, a value 
of x = 25% is chosen to ensure minimal anisotropy. 


4.3. Discrete domain refining 

The three criteria that drive the domain homogene- 
ity and isotropy are the volume fraction, the cardinal 
number and the mean square difference of the contact 
orientation packet. 

This section deals with two questions: 


e Do the three criteria converge if the number of dis- 
crete element per unit of volume increases ? 


(a) k = 0% 


(b) k = 25% 


Figure 6: 3D histograms of the orientation set 
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Figure 7: Evolution of the mean square difference e pa- 
rameter of the sampling distribution of the contact ori- 
entation packet versus radius dispersion x 


e Inthis case, is it possible to define the right number 
of discrete elements that allows the simulation of 
an homogeneous and isotropic medium ? 


First, the meaning of refining must be clarified. Refining 
consists of increasing the number of discrete elements 
per unit of volume. To study the influence of refining 
on the three criteria, 22 packing domains of identical 
bounding volume* were built. To permit a statistical 
processing, 5 packing domains are built for each dis- 
crete element number value. Therefore, a total of 110 
packing domains were analyzed. 

Figures 8a, 8b and 8c show the influence of refining 
on the volume fraction, cardinal number and isotropy. 
These values are extracted from the discrete domains 
built with a dispersion radius of k = 25% (correspond- 
ing to the conclusion of the previous section ). Small 
differences could be accepted. The criteria could be 
classified by order of importance : 


1. Isotropy is considered as the most important. This 
criterion highly influences the discrete sample me- 
chanical behavior. 

2. The cardinal number and the volumic fraction are 
less important. They allow to check the validity 
of the compact domain. Small differences between 
the sample values and the RCP values established 
in the literature are accepted. 


The cardinal number converges to a limit value close to 
6.2. For Volume fraction and isotropy slight variations 
are still observed for high number of discrete elements. 

The discrete element number increases with refining. 
A high refining level brings down the computational 
performances. Therefore, a compromise must be made 
between performance and precision. For the next sec- 
tion, anumber of 10 000 discrete elements is considered 
as sufficient for an acceptable level of precision. For this 
value, the geometrical anisotropy criterion is lower than 
0.0032, the coordination number is higher than 6.2 and 
the volume fraction around 0.635. 10 000 discrete el- 
ements in a 3D square domain gives the ratio between 
the micro-scale and the macro-scale to obtain isotropy 
and homogeneity at the macroscale; this ratio is around 
10000!/ x 21.5. 


4.4. Conclusion 

The method presented in this section shows how to 
verify the initial compacted domain against the geomet- 
rical criteria. The geometrical criteria do not depend on 


4Bounding volume is defined as the volume of the box of minimal 
size containing the discrete domain. 


the discrete element mechanical interaction laws. The 
presented method can be applied to various discrete el- 
ement models and applications. 

An original method has been presented to quantify 
the isotropy of a discrete element set. This method is 
based on the classification of bond directions by using a 
platonic solid (geode). 

Finally, a ratio of "macroscale over microscale" about 
20 has been chosen, from the convergence curves, that 
control the geometrical properties of the discrete ele- 
ment set. 


5. Elastic calibration 


The previous section introduced a methodology to 
build an initial compact discrete domain that ensures an 
adequate, homogeneous, isotropic and geometric orga- 
nization. Then, cohesive beams (see section 2) are intro- 
duced by creating a beam at each contact between two 
discrete elements. In further simulations, the contacts 
are not taken in account; only cohesive beam interac- 
tions between discrete elements are considered. 

To study the influence of cohesive beam bond pa- 
rameters on macroscopic elastic behavior, a parametric 
study using numerical quasi-static uni-axial tensile test 
is used. 


5.1. Quasi-static tensile test description 


A perfectly homogeneous, isotropic, elastic material 
is characterized by the Young’s modulus and the Pois- 
son’s ratio. For real materials, these parameters are gen- 
erally determined by quasi-static tensile tests. These ex- 
perimental procedures can be also applied to a numeri- 
cal sample. 


5.1.1. From discrete to continuous geometry 

To compute the macroscopic Young’s modulus and 
Poisson’s ratio, a perfect 3D continuous geometry is as- 
sociated to the compact discrete domain. This perfect 
geometry is the bounding shape of the compacted dis- 
crete domain. The discrete elements belonging to the 
domain boundaries are marked to compute the perfect 
geometry dimensions. With the cylinder shape, three 
discrete element set are marked (see figure 9): 


e the "xMax" and "xMin" set are associated to faces 
with normal X and —X. 


e the "radius" set is associated with the cylinder cir- 
cumference. 
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Figure 8: Evolution of 3 geometrical criteria versus refining (for radius dispersion x = 25%) 


The Perfect cylinder dimensions are computed as : 


Nxmax Nxmin 
= —= 1 — | = 
Lm = 2R + OG; + OG; X 
xMax =] Nxmin i=0 
(12) 
1 Nradius = = 
Ru = R + — > (OG;.Y)2 + (OG;.Z)2 (13) 
Nadius i=l 
where: 


e Ly and Ry are the perfect cylinder length and ra- 
dius. 


© NMax, Nemin and Nadius are the number of dis- 
crete elements belonging to "xMax", "xMin" and 
"radius" sets. 


— 
e OG; is the position of the discrete element gravity 
center. 


e R is the average discrete element radius over the 
entire domain. 


5.1.2. Loading 


To ensure a quasi-static tensile test, the loading force 
acting on the discrete element set "xMin" and "xMax", 
are progressively applied (linear ramp) and stabilized. 
The sum of forces acting on "xMax" and "xMin" are 


— —_ | 
denoted by Fmax and Fyyin. These two forces, acting 


=> a = — 

along X for Fyyax and —X for Fxmin, have equal norms 
and are opposed. To check the quasi-static properties, 
the kinetic and deformation energies are computed and 
stored during the numerical test. Figure 10 confirms that 
a progressive loading gives a negligible kinetic energy 
and ensures a quasi-static aspect of the simulation. 


Figure 9: Perfect cylinder associated with a discrete do- 
main with x = 25% 
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Figure 10: Kinetic and deformation energy during a 
quasi-static tensile test (computed with a time step At = 
3.1077s and a number of iteration it = 100 000) 


5.1.3. Young’s modulus and Poisson’s ratio computa- 
tion 
The Young’s modulus and the Poisson’s ratio can be 
easily determined for the cylinder sample by using the 
material strength analytical formulations: 


F/S m, 
PM = Kou Lu, 
ARu/Rm, 
gee sacs St 1 
YM TANT (15) 
(16) 


where: 


e Lm, Rm, and S m, are the initial bounding cylinder 
dimensions (respectively: length, radius, and sec- 
tion). 


e Ey and vy are the macroscopic Young modulus 
and Poisson ratio. 


e F is the normal force. 


The explicit numeric schemes are not well-adapted to 
the quasi-static simulation. The system vibrates around 
the static solution. To allow a convergence, a pure nu- 
merical damping factor is introduced in the numerical 
scheme as described by Tchamwa and Wielgosz (38). 
This is a decentered explicit integration scheme that al- 
low high frequency dissipation. This scheme is very 
similar to the Velocity Verlet algorithm. The dissipation 
is controlled with a single parameter y. Only the second 
time derivative equality is modified. The equation 6 and 
9 become : 


> > d > > 
B+ dr) = B+ eS (BD + BE+dy) (M) 


d 
gerdd = an +65 GOAD) (18) 


The value 6 = 1.3 is used to allow a high convergence 
rate to the static solution. 


5.2. Parametric study 


The cohesive beam bond is defined by four parame- 
ters: 


e two geometrical parameters: length L, and radius 


Ty. 


e two mechanical parameters: Young’s modulus E, 
and Poisson’s ratio v,. 
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The value of the cohesive beam bond length depends on 
the distance between discrete element centers and is not 
a free parameter. The three others parameters are free 
and must be quantified. 

The adimensional cohesive beam radius parameter 
noted À, will be preferred to the beam radius. It al- 
lows a definition that does not depend on the discrete 
element sizes. It is defined as the ratio between the co- 
hesive beam radius and the average discrete element ra- 
dius R. Note that this value is the same for all the co- 
hesive beams; consequently, all the cohesive beam radii 
are equal. 


5.2.1. Microscopic Poisson’s ratio influence 

Figure 11 shows the evolution of the macroscopic 
Young’s modulus E and Poisson’s ratio vy for the dif- 
ferent microscopic Poisson’s ratio v, values in the range 
[0, 1/2]. It is observed that the microscopic Poisson’s 
ration y, does not influence the macroscopic Young’s 
modulus Ey and the macroscopic Poisson’s ratio vm 
significantly. 

Equations 3 and 4 show that the microscopic shear 
modulus G,, and consequently the microscopic Pois- 
son’s ratio v, plays a role in local torsion loading (in 
a cohesive beam). During a quasi-static tensile test, the 
contribution of the local torsion energy is minor (see fig- 
ure 12). Each energy represented on this figure are the 
sum of the local elastic energy stored by the cohesive 
beams. The total elastic energy is split up into : 


Tension energy characterized by the sum of the cohe- 
sive beam works of the normal forces. 


Bending energy characterized by the sum of the cohe- 
sive beam works of the bending torques. 


Torsion energy characterized by the sum of the cohe- 
sive beam works of the torsion torques. 


To summarize, the influence of the microscopic Pois- 
son’s ratio is negligible. Consequently, the microscopic 
Poisson’s ratio values can be chosen arbitrarily. For the 
rest of the study a value of 0.3 is chosen. 


5.2.2. Microscopic Young’s modulus influence 

Figure 13 shows the evolution of the macroscopic pa- 
rameters as a function of microscopic Young modulus 
for different values of the microscopic radius ratio. The 
next table gives an outline of these evolutions. 


Macroscopic Macroscopic Young’s Macroscopic Poisson’s 
parameters modulus Ey ratio vy 
Functions f 2 
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Figure 13: Microscopic Young’s modulus £, influence on the macroscopic parameters Ey and vy 
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Figure 11: Influence of v, on Ey and vm 
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Figure 12: Energy breakdown of total elastic energy 
stored by cohesive beams for the quasi-static tensile test 
(computed with a time step At = 3.10~'s and a number 
of iteration it = 100 000) 
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5.2.3. Microscopic radius ratio influence 

Figure 14 shows the evolution of the macroscopic pa- 
rameters as a function of the microscopic radius ratio ^, 
for different values of the microscopic Young’s modu- 
lus. The next table gives an outline of these evolutions. 


Macroscopic Macroscopic Young’s Macroscopic Poisson’s 
parameters modulus Ey ratio Yy 


Functions Vu = 
decreasing quadratic 
function 
14b 


increasing quadratic 
function. 
14a 


Evolution 


Figures 


5.3. Calibration method 


Section 5.2 has described the influences of the mi- 
croscopic parameters v,, E, and 7, on the macroscopic 
parameters Em and vy. These influences will be used 
to develop a calibration methodology in two steps: the 
calibration of microscopic radius ratio, then the calibra- 
tion of microscopic Young’s modulus. Section 5.2.1 has 
shown that the influence of the microscopic Poisson’s 
ratio v, is negligible. Its value is arbitrarily fixed at 0.3. 


5.3.1. Microscopic radius ratio F, calibration 

Section 5.2.2 have shown that the influence of the mi- 
croscopic Young’s modulus E, is very small on macro- 
scopic Poisson ratio vy (figure 13b). 

The first calibration step considers that the macro- 
scopic Poisson’s ratio vy does not depend on the mi- 
croscopic Young’s modulus E,,. A single macroscopic 
Poisson’s ratio vy value is associated with each micro- 
scopic radius ratio f,. The figure 15 shows the evolution 
of the macroscopic Poisson’s ratio vy as a function of 
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Figure 14: Microscopic radius ratio ^, influence on the macroscopic parameters Ey and vm 
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Figure 15: Radius ratio 7, calibration for fused silica 


the microscopic radius ratio ,. This description allows 
the deduction of a microscopic radius ratio Ê, value cor- 
responding to a desired-value of the macroscopic Pois- 
son’s ratio vy. Figure 15 shows an application for the 
silica glass material. In this case, a microscopic radius 
ratio value off," ~ 0.71 is found to correspond to 
the silica Poisson’s ratio value ye z 0.17. 


5.3.2. Microscopic Young modulus E, calibration 
Figure 16 shows the evolution of the macroscopic 
Young’s modulus Ey as a function of the microscopic 
Young’s modulus £, for a radius ratio (evaluated in the 
previous section) of ,°""“¢ = 0.71. This evolution al- 
lows the deduction of a microscopic Young’s modulus 
value for a desired-value of the macroscopic Young’s 
modulus. In the case of silica glass, the microscopic 


Young’s modulus Es" ~ 265 GPa is found to corre- 
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Figure 16: Microscopic Young modulus E, calibration 
of fused silica 


spond to the silica glass Young’s modulus value E ae x 
72.5 GPa (see figure 16). 


5.4. Study of assembly dependency 


To apply the micro beam model to any material ge- 
ometry, it must be verified that the calibration results 
do not depend on the number of discrete elements in a 
given material volume. To check this property, many 
discrete samples (with similar bounding dimensions) 
were built with an increasing number of discrete ele- 
ment (see figures 18). The samples satisfy the criteria 
established in section 4. To take into account the vari- 
ability of the sample geometry, four different samples 
were built with the same discrete element number. The 
figure 17 shows the evolution of the macroscopic pa- 
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Figure 17: Evolution of the macroscopic parameters Em 


and vy as a function of the discrete element number. 


rameters Ey and vy as a function of the discrete ele- 
ment number. 

It appears that for a number of discrete element over 
7 500, the macroscopic Young modulus Ey fluctuates 
around 7% and the macroscopic Poisson ratio around 
YM 5%. 

Hentz (7) shown that the Liao calibration methodol- 
ogy for the dual spring model (39) gives a dispersion 
around 28% for the Young modulus and 16% for the 
Poisson’s ratio. To improve the accuracy, Hentz has 
introduced an energy criterion to reduce the dispersion 
around 12% and 10%. But this criterion is assembly 
dependent and must be computed for each sample. 

The beam cohesive model associated to the com- 
paction criteria presented in section 4 allow a better pre- 
cision without any re-computation. 


5.5. Validation tests 


The previous subsections show a methodology to cal- 
ibrate the microscopic parameters from the macroscopic 
elastic parameters values. Tables 1 and 2 summarize the 
results obtained for the silica glass material. 

These microscopic elastic parameters are used to 
build a cylindrical numeric sample of silica glass. This 
sample is submitted to quasi-static tensile, bending and 
torsion testing, in which the "xMin" discrete element set 
is fixed and "xMax" set is loaded. To reproduce a quasi- 
static aspect, the loads are applied gradually (cf section 
5.1.2). 

Free face displacement and rotation given by the nu- 
merical simulations are compared to the results given 
by the strength of material theory. Table 3 summarizes 
the differences as a percent between numerical and the- 
oretical results. The higher difference is obtained for 
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Young’s Modulus 
Eu = 123 GPa 


Table 1: Macroscopic silica glass elastic values 


Young's Modulus 
E, = 265 GPa 


Table 2: Microscopic silica glass elastic values 


the torsion test and is less than seven percent. These er- 
rors are considered acceptable. The causes are multiple 
and are inherent to the random discrete element posi- 
tioning: imperfect loading and imprecise measurement 
of boundary geometry. 


6. Dynamic calibration 


The previous section deals with the elastic calibra- 
tion. This method allows to the calibration of the three 
microscopic elastic parameters (Young’s modulus, ra- 
dius and Poisson’s ratio) to obtain the elastic behavior 
at the macroscopic scale. To quantitatively simulate dy- 
namic phenomena likes cracks or impacts, it is also nec- 
essary to calibrate the microscopic mass parameters. 

The discrete elements mass parameters (inertia ma- 
trix and masses) depend on the discrete element volume 
and density. The discrete element geometries are ini- 
tially determined. Therefore, only the discrete element 
density can be adjusted. In the same way as the elas- 
tic parameters, the microscopic density can be differ- 
ent from the macroscopic density to correct the voids 
between the discrete elements in a compacted domain. 
Therefore, the density will be considered as the only 
free parameter. 

A very simple calibration criterion is chosen. This 
criterion ensures mass equality between the discrete and 
continuous domains: 


(19) 


where: 


e p, and V, are the discrete element density and vol- 
ume. 


e pm and Vy are the continuous density and volume. 
The continuous domain dimensions are computed 
as presented in section 5.1.1. 


(a) 200 
discrete element 


(b) 2 000 
discrete element 


(c) 20 000 
discrete element 


Figure 18: Snapshot of discrete samples with increasing fineness 


Tse ting er 
| Criteria || Free face normal displacement | Free face tangential deflection 
| Difference || 1.20 % 4.16 % 6.13 % 


Table 3: Comparison of the numerical and theoretical results for the quasi-static tensile, bending and torsion tests 
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Figure 19: Energy breakdown of total elastic energy 
stored by cohesive beams for a complex dynamic test 
(computed with a time step At = 3.10~'s and a number 
of iteration it = 50 000) 


However, this method does not ensure inertia equality. 
This difference is supposed as negligible. This assump- 
tion is based on: 


1. a classical strength material hypothesis, i.e., for 
beam transverse vibrations, the cross section rota- 
tion energies are negligible compared to the trans- 
lation energies. 

2. the sum of the local torsion energies is negligible 
in the quasi-static (see figure 12) and dynamic (see 
figure 19) DEM tests. 


This section deals with the validation of this hypothe- 
sis by verifying the calibration method by dynamic ten- 
sile, bending, torsion and impact tests on the discrete 
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domain (cf figure 9). The numerical results are com- 
pared to the analytical results given for the associated 
continuous domain. 

The macroscopic characteristics correspond to silica 
glass, with a Young’s modulus, Poisson’s ratio, and 
density equal to Ey = 72.5 GPa, vy = 0.17 and 
pm = 2201 kg/m°. The microscopic elastic parameters 
correspond to those deduced in section 5. The micro- 
scopic density is deduced from the calibration method 
based on mass equality between the discrete and the 
continuous domains. The continuous domain dimen- 
sions are computed by the method presented in section 
5.1.1. 

For the dynamic tensile, bending and torsion tests, 
the boundary conditions and loadings are : "xMin" fixed 
and "xMax" loaded. The force acting on "xMax" is ap- 
plied gradually and then suddenly set to 0. This loading 
allows an excitation of the transient response of the dis- 
crete domain. 

The oscillation periods are computed from a frequen- 
cial analysis (FFT) of the average free face position and 
angular orientation. The numerical results are compared 
to the theoretical results given by the vibration of the 
continuous system analytical models (40, §4). 

For the impact tests an initial velocity on X is applied 
on the "xMax" face to generate a mechanical shock 
wave. The "xMin" average velocity on X is measured at 
each time step. So, it is possible to observe the moment 
when the mechanical wave reaches the "xMin" face (see 
figure 20). The mechanical wave celerity is deduced 
from the elapsed time corresponding to the mechanical 
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Figure 20: "xMin" and "xMax" face average velocities 
as a function of time 
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Figure 22: Cross section computational improvement 


wave flight between the "xMin" and "xMax" faces (fig- 
ures 21). 


Table 4 shows the comparisons between the simula- 
tions and theoretical results. Tensile, torsion and impact 
tests show a very good agreement with the analytical re- 
sults. The bending test shows less precision. In fact, the 
error depends on the cross-sectional moment of inertia 
computation. In addition, an error in the measurement 
of the radius has a high influence on results. Figure 22 
shows the "xMax" discrete sample face. Following the 
definition established in section 5.1.1, the sample sec- 
tion is the "Max section". However, the right section, 
in which internal forces exist, is closer to the "Effective 
Section". With this new definition of the sample radius, 
the error between the numerical and the theoretical re- 
sults for dynamic bending tests is around 1.8 %. At the 
moment, this new definition is too experimental. Its im- 
pact on the elastic calibration and the dynamic method 
needs to be further explored. 
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7. Conclusion 


A methodology to verify the initial compact domain 
using geometrical criteria : the isotropy, the cardinal 
number, the volume fraction and the fineness has been 
presented. The most important criterion that further im- 
pacts mechanical properties is the geometrical isotropy. 
An original method, based on geometrical and statistical 
analysis, is presented to ensure a good level of isotropy. 
A simple fineness criterion is also presented to ensure a 
stable geometrical criterion. 

A methodology has been presented to obtain the de- 
sired value of the Young’s modulus and the Poisson’s 
ratio at the macroscopic scale for a 3D spherical discrete 
element bonded by microscopic beams. This calibration 
is based on a parametric study of the microscopic and 
macroscopic elastic parameters. The deduced curves 
can be used as an "abacus" to calibrate the elastic pa- 
rameters easily. 

A calibration method for the dynamic parameters has 
been presented. The discrete element density is com- 
puted to ensure equality between the continuous and the 
total discrete domain mass. 

The calibration results do not depend on the discrete 
element size. This important property validates the in- 
terest of the micro beam model and the proposed cali- 
bration method. 

Numerical samples of silica glass has been calibrated. 
This sample has been tested under tensile, bending, and 
torsion quasi-static and dynamic tests. The results show 
good agreement with the strength of material theory. 

In conclusion, with the proposed methodology, a dis- 
crete element model for homogeneous and isotropic ma- 
terials is obtained with good quantitative results. DEM 
has often been used as a qualitative tool to understand 
complex phenomena such as wear, fracturing or impact. 
This work is a first step to propose a quantitative numer- 
ical tool that will be able to propose predictive models 
for these classes of problems that have no predictive nu- 
merical model presently. 


8. Acknowledgements 


This work is supported by the Conseil Régional 
d’ Aquitaine and is performed in the framework of the 
Etude et Formation en Surfacage Optique (EFESO) 
project. 

The developments carried out during this project 
have been implemented in the Granoo Project, see 
http://www.granoo.org for details. 


(a) t x 3.8 us (b) t x 7.6 us 


(c) t x 11.4 us 


(d) t x 15.2 us (e) t x 19 us 


Figure 21: Snapshots of mechanical wave propagation 
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Table 4: Comparison of the numerical and theoretical results for the dynamic tensile, bending, torsion and impact 
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